ARGUS distribution (argus)#
Goal: build intuition, derive key formulas, and implement simulation/visualization for the continuous ARGUS distribution as implemented in scipy.stats.argus.
The ARGUS distribution is most famous in particle physics as a parametric model for background shapes near a hard kinematic endpoint.
import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
import scipy
from scipy import special
from scipy.optimize import minimize_scalar
from scipy.stats import argus, chi2, norm
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
rng = np.random.default_rng(42)
print("numpy:", np.__version__)
print("scipy:", scipy.__version__)
numpy: 1.26.2
scipy: 1.15.0
1) Title & classification#
Name:
argus(ARGUS distribution)Type: continuous
Support (standardized): (x \in (0, 1))
The density is 0 at the endpoints, so people often write ([0,1]) informally.
Parameter space: shape parameter (\chi > 0)
SciPy also supports
locandscale, shifting/scaling the support to ((\text{loc},, \text{loc}+\text{scale})).
2) Intuition & motivation#
What it models#
A typical use case is a background distribution for a quantity with a hard upper endpoint (e.g., an invariant mass that cannot exceed a known limit).
The factor (\sqrt{1-x^2}) creates a phase-space-like suppression as (x \to 1).
The exponential factor (\exp{-\tfrac{\chi^2}{2}(1-x^2)}) controls how strongly mass piles up near the endpoint.
Real-world use cases#
High energy physics: background modeling near kinematic limits (the classic “ARGUS function”).
Endpoint distributions: any normalized measurement bounded above, where the density falls to 0 at the maximum.
Relations to other distributions#
A very useful transformation connects ARGUS to a truncated Gamma distribution:
[ Y ;:=; \frac{\chi^2}{2}(1-X^2) \in (0, \chi^2/2). ]
Then (Y) has a density proportional to (\sqrt{y} e^{-y}), i.e. a (\text{Gamma}(3/2,,1)) conditioned on ([0,\chi^2/2]).
Limiting behavior:
As (\chi \downarrow 0): the exponential term (\to 1), and the density approaches (f(x) \propto x\sqrt{1-x^2}).
As (\chi \to \infty): the distribution concentrates near (x \approx 1).
3) Formal definition#
PDF#
For (0<x<1) and (\chi>0), the standardized ARGUS density is
[ f(x;\chi) = \frac{\chi^3}{\sqrt{2\pi},\Psi(\chi)}, x,\sqrt{1-x^2},\exp\left{-\frac{\chi^2}{2}(1-x^2)\right}, ]
where
[ \Psi(\chi) = \Phi(\chi) - \chi,\phi(\chi) - \tfrac{1}{2}, ]
and (\Phi) and (\phi) are the CDF and PDF of (\mathcal N(0,1)).
A numerically stable equivalent form is
[ \Psi(\chi) = \tfrac{1}{2},P\left(\tfrac{3}{2}, \tfrac{\chi^2}{2}\right), ]
where (P(a,x)=\gamma(a,x)/\Gamma(a)) is the regularized lower incomplete gamma function.
CDF#
SciPy’s implementation uses a particularly clean survival function:
[ \bar F(x;\chi) = 1 - F(x;\chi) = \frac{\Psi\big(\chi\sqrt{1-x^2}\big)}{\Psi(\chi)}, ]
so
[ F(x;\chi) = 1 - \frac{\Psi\big(\chi\sqrt{1-x^2}\big)}{\Psi(\chi)}. ]
Location/scale#
If (Z\sim\text{ARGUS}(\chi)) on ((0,1)), then (X=\text{loc}+\text{scale}\cdot Z) has support ((\text{loc},, \text{loc}+\text{scale})).
def argus_Psi(chi: np.ndarray | float) -> np.ndarray:
"""Psi(chi) used by the ARGUS distribution.
We use the regularized incomplete gamma form for numerical stability:
Psi(chi) = 0.5 * P(3/2, chi^2/2).
"""
chi = np.asarray(chi, dtype=float)
return 0.5 * special.gammainc(1.5, chi**2 / 2)
def argus_pdf(x: np.ndarray | float, chi: float) -> np.ndarray:
x = np.asarray(x, dtype=float)
chi = float(chi)
if chi <= 0:
return np.full_like(x, np.nan, dtype=float)
Psi = argus_Psi(chi)
norm_const = chi**3 / (np.sqrt(2 * np.pi) * Psi)
y = 1.0 - x * x
base = norm_const * x * np.sqrt(np.clip(y, 0.0, None)) * np.exp(-0.5 * chi**2 * y)
return np.where((x > 0) & (x < 1), base, 0.0)
def argus_cdf(x: np.ndarray | float, chi: float) -> np.ndarray:
x = np.asarray(x, dtype=float)
chi = float(chi)
if chi <= 0:
return np.full_like(x, np.nan, dtype=float)
Psi_chi = argus_Psi(chi)
t = chi * np.sqrt(np.clip(1.0 - x * x, 0.0, None))
sf = argus_Psi(t) / Psi_chi
cdf = 1.0 - sf
cdf = np.where(x <= 0, 0.0, cdf)
cdf = np.where(x >= 1, 1.0, cdf)
return cdf
# Quick consistency check vs SciPy
xgrid = np.linspace(0, 1, 400)
chi_test = 2.5
max_pdf_err = np.max(np.abs(argus_pdf(xgrid, chi_test) - argus.pdf(xgrid, chi_test)))
max_cdf_err = np.max(np.abs(argus_cdf(xgrid, chi_test) - argus.cdf(xgrid, chi_test)))
print("max |pdf - scipy|:", max_pdf_err)
print("max |cdf - scipy|:", max_cdf_err)
max |pdf - scipy|: 1.7763568394002505e-15
max |cdf - scipy|: 3.552713678800501e-15
4) Moments & properties#
Mean and variance#
SciPy uses (and we can derive) the following closed forms.
Let (I_1) be the modified Bessel function of the first kind. Define (z=\chi^2/4).
[ \mathbb E[X] = \sqrt{\frac{\pi}{8}},\frac{\chi,e^{-z} I_1(z)}{\Psi(\chi)}. ]
For the second moment: [ \mathbb E[X^2] = 1 - \frac{3}{\chi^2} + \frac{\chi,\phi(\chi)}{\Psi(\chi)}. ]
Then (\operatorname{Var}(X)=\mathbb E[X^2]-\mathbb E[X]^2).
Skewness and kurtosis#
Skewness and (excess) kurtosis are defined via central moments: [ \gamma_1 = \frac{\mu_3}{\sigma^3},\qquad \gamma_2 = \frac{\mu_4}{\sigma^4} - 3. ]
For ARGUS, these don’t simplify nicely to a short expression; they’re typically computed numerically.
MGF / characteristic function#
Because (X\in(0,1)) is bounded, both exist for all real arguments: [ M_X(t)=\mathbb E[e^{tX}],\qquad \varphi_X(\omega)=\mathbb E[e^{i\omega X}]. ]
There is no widely used simple closed form; numerical quadrature (or a moment series) is standard.
Entropy#
The differential entropy is [ H(X) = -\int_0^1 f(x;\chi),\log f(x;\chi),dx. ]
SciPy provides argus.entropy(chi).
def argus_mean(chi: np.ndarray | float) -> np.ndarray:
chi = np.asarray(chi, dtype=float)
Psi = argus_Psi(chi)
z = chi**2 / 4
return np.sqrt(np.pi / 8) * chi * special.ive(1, z) / Psi
def argus_E_x2(chi: np.ndarray | float) -> np.ndarray:
"""E[X^2] with a small-chi branch for numerical stability (mirrors SciPy)."""
chi = np.asarray(chi, dtype=float)
out = np.empty_like(chi)
mask = chi > 0.1
if np.any(mask):
c = chi[mask]
out[mask] = 1 - 3 / c**2 + c * norm.pdf(c) / argus_Psi(c)
if np.any(~mask):
c = chi[~mask]
# series approximation for small chi (from SciPy's implementation)
coef = [-358 / 65690625, 0, -94 / 1010625, 0, 2 / 2625, 0, 6 / 175, 0, 0.4]
out[~mask] = np.polyval(coef, c)
return out
def argus_var(chi: np.ndarray | float) -> np.ndarray:
m = argus_mean(chi)
return argus_E_x2(chi) - m**2
chis = np.array([0.2, 0.5, 1.0, 2.5, 6.0])
m_formula = argus_mean(chis)
v_formula = argus_var(chis)
m_scipy, v_scipy, s_scipy, k_scipy = argus.stats(chis, moments="mvsk")
print("chi mean(formula) mean(scipy) var(formula) var(scipy) skew kurt(excess)")
for chi, mf, ms, vf, vs, ss, ks in zip(chis, m_formula, m_scipy, v_formula, v_scipy, s_scipy, k_scipy):
print(f"{chi:4.1f} {mf:12.6f} {ms:11.6f} {vf:12.6f} {vs:11.6f} {ss:7.4f} {ks:11.4f}")
# MGF/CF: demonstrate numerical quadrature vs Monte Carlo
chi0 = 2.5
x_mc = argus.rvs(chi0, size=150_000, random_state=rng)
ts = np.array([-3.0, 0.0, 3.0])
ws = np.array([0.0, 10.0, 20.0])
print("\nMGF M(t)=E[e^{tX}] at a few t:")
for t in ts:
mc = np.mean(np.exp(t * x_mc))
quad = argus.expect(lambda x, t=t: np.exp(t * x), args=(chi0,))
print(f"t={t:>5.1f} MC={mc:.6f} quad={quad:.6f}")
print("\nCharacteristic function φ(ω)=E[e^{iωX}] at a few ω:")
for w in ws:
mc = np.mean(np.exp(1j * w * x_mc))
quad = argus.expect(lambda x, w=w: np.exp(1j * w * x), args=(chi0,))
print(f"w={w:>5.1f} MC={mc:.6f} quad={quad:.6f}")
print("\nEntropy at chi=2.5:", argus.entropy(chi0))
chi mean(formula) mean(scipy) var(formula) var(scipy) skew kurt(excess)
0.2 0.590227 0.590227 0.053005 0.053005 -0.2966 -0.7982
0.5 0.596428 0.596428 0.052891 0.052891 -0.3230 -0.7800
1.0 0.618703 0.618703 0.052157 0.052157 -0.4204 -0.6939
2.5 0.762779 0.762779 0.035555 0.035555 -1.1882 1.0606
6.0 0.956717 0.956717 0.001359 0.001359 -1.8768 5.8944
MGF M(t)=E[e^{tX}] at a few t:
t= -3.0 MC=0.123343 quad=0.123745
t= 0.0 MC=1.000000 quad=1.000000
t= 3.0 MC=11.224313 quad=11.216045
Characteristic function φ(ω)=E[e^{iωX}] at a few ω:
w= 0.0 MC=1.000000+0.000000j quad=1.000000
w= 10.0 MC=-0.252645+0.292744j quad=-0.253406
w= 20.0 MC=0.133822-0.107030j quad=0.133083
Entropy at chi=2.5: -0.48713485149744684
/home/tempa/miniconda3/lib/python3.12/site-packages/scipy/integrate/_quadpack_py.py:606: ComplexWarning:
Casting complex values to real discards the imaginary part
5) Parameter interpretation (shape changes)#
The single shape parameter (\chi) controls how quickly the density falls away from the endpoint.
Small (\chi): the exponential term is weak, and the shape is dominated by (x\sqrt{1-x^2}).
Large (\chi): the exponential term strongly favors (x) near 1, creating a sharp peak close to the endpoint.
A convenient closed-form for the mode (maximum of the pdf) comes from differentiating (\log f(x)):
[ \text{mode}(X)^2 = \frac{\chi^2 - 2 + \sqrt{\chi^4 + 4}}{2\chi^2}. ]
def argus_mode(chi: float) -> float:
chi = float(chi)
if chi <= 0:
return np.nan
u = (chi**2 - 2 + np.sqrt(chi**4 + 4)) / (2 * chi**2)
return float(np.sqrt(u))
for chi in [0.2, 1.0, 2.5, 6.0]:
print(f"chi={chi:>4}: mode≈{argus_mode(chi):.4f}")
chi= 0.2: mode≈0.7106
chi= 1.0: mode≈0.7862
chi= 2.5: mode≈0.9300
chi= 6.0: mode≈0.9864
6) Derivations#
6.1 Derivation of the mean (why Bessel functions appear)#
Start from ( \mathbb E[X] = \int_0^1 x,f(x;\chi),dx )
so the integral of interest is
[ I = \int_0^1 x^2\sqrt{1-x^2},\exp\left{-\frac{\chi^2}{2}(1-x^2)\right},dx. ]
Use the substitution (x=\cos\theta) with (\theta\in[0,\pi/2]). Then (1-x^2=\sin^2\theta), (dx=-\sin\theta,d\theta), and
[ I = \int_0^{\pi/2} \cos^2\theta,\sin^2\theta,\exp\left{-\frac{\chi^2}{2}\sin^2\theta\right}d\theta. ]
Rewrite (\sin^2\theta = \tfrac{1-\cos 2\theta}{2}) and set (\varphi=2\theta), which yields an integral of the form (\int_0^{\pi} e^{z\cos\varphi}\cos(n\varphi)d\varphi), whose value is (\pi I_n(z)).
After simplification you arrive at
[ I = \frac{\pi}{2\chi^2} e^{-\chi^2/4} I_1(\chi^2/4). ]
Multiplying by the normalization constant (\chi^3/(\sqrt{2\pi}\Psi(\chi))) gives
[ \mathbb E[X] = \sqrt{\frac{\pi}{8}},\frac{\chi,e^{-\chi^2/4}I_1(\chi^2/4)}{\Psi(\chi)}. ]
6.2 Derivation of (\mathbb E[X^2]) via the truncated-Gamma representation#
Let (Y=\tfrac{\chi^2}{2}(1-X^2)). Then (Y) is Gamma((3/2,1)) conditioned on ([0,\chi^2/2]).
Since (X^2=1-2Y/\chi^2):
[ \mathbb E[X^2] = 1 - \frac{2}{\chi^2},\mathbb E[Y\mid Y\le \chi^2/2]. ]
For a Gamma((\alpha,1)) truncated at (b), (\mathbb E[Y\mid Y\le b] = \gamma(\alpha+1,b)/\gamma(\alpha,b)).
Using the incomplete-gamma recurrence simplifies the ratio and leads to
[ \mathbb E[X^2] = 1 - \frac{3}{\chi^2} + \frac{\chi,\phi(\chi)}{\Psi(\chi)}. ]
6.3 Likelihood for i.i.d. data#
For data (x_1,\dots,x_n\in(0,1)), the log-likelihood (standardized) is
[ \ell(\chi) = n\left(3\log\chi - \log\Psi(\chi) - \tfrac{1}{2}\log(2\pi)\right)
\sum_{i=1}^n\left(\log x_i + \tfrac{1}{2}\log(1-x_i^2) - \tfrac{\chi^2}{2}(1-x_i^2)\right). ]
A handy derivative identity is (\Psi’(\chi)=\chi^2\phi(\chi)), so you can optimize (\ell(\chi)) with gradient-based methods.
def argus_loglik(chi: float, x: np.ndarray) -> float:
"""Log-likelihood for standardized ARGUS(chi) given x in (0,1)."""
chi = float(chi)
if chi <= 0:
return -np.inf
x = np.asarray(x, dtype=float)
if np.any((x <= 0) | (x >= 1)):
return -np.inf
Psi = argus_Psi(chi)
y = 1.0 - x * x
ll = (
x.size * (3 * np.log(chi) - 0.5 * np.log(2 * np.pi) - np.log(Psi))
+ np.sum(np.log(x) + 0.5 * np.log1p(-x * x) - 0.5 * chi**2 * y)
)
return float(ll)
def argus_loglik_grad(chi: float, x: np.ndarray) -> float:
"""Gradient of the log-likelihood.
Uses Psi'(chi) = chi^2 * phi(chi).
"""
chi = float(chi)
if chi <= 0:
return np.nan
x = np.asarray(x, dtype=float)
if np.any((x <= 0) | (x >= 1)):
return np.nan
n = x.size
Psi = argus_Psi(chi)
Psi_prime = chi**2 * norm.pdf(chi)
y = 1.0 - x * x
return float(n * (3 / chi - Psi_prime / Psi) - chi * np.sum(y))
# Simulate data and compute the MLE
chi_true = 2.5
x = argus.rvs(chi_true, size=800, random_state=rng)
obj = lambda c: -argus_loglik(c, x)
res = minimize_scalar(obj, bounds=(1e-6, 30.0), method="bounded")
chi_hat = float(res.x)
print("chi_true:", chi_true)
print("chi_hat (MLE):", chi_hat)
print("grad at chi_hat:", argus_loglik_grad(chi_hat, x))
# Likelihood ratio test: H0: chi = chi0 vs H1: chi free
chi0 = 1.0
lr_stat = 2 * (argus_loglik(chi_hat, x) - argus_loglik(chi0, x))
p_value = chi2.sf(lr_stat, df=1)
print("\nLRT statistic:", lr_stat)
print("approx p-value (chi-square df=1):", p_value)
chi_true: 2.5
chi_hat (MLE): 2.5179774750503756
grad at chi_hat: -2.3556725523121713e-05
LRT statistic: 386.62155391177384
approx p-value (chi-square df=1): 4.5017064542421926e-86
7) Sampling & simulation (NumPy-only)#
Below is a NumPy-only sampler for standardized ARGUS((\chi)). It mirrors SciPy’s piecewise strategy:
Small (\chi): rejection sampling using the (\chi\to 0) base density (g(x)=3x\sqrt{1-x^2}).
Moderate (\chi): rejection sampling with a proposal density (g(x)\propto x\exp{-\tfrac{\chi^2}{2}(1-x^2)}).
Large (\chi): use the truncated Gamma representation of (Y=\tfrac{\chi^2}{2}(1-X^2)).
This is not the only way to sample ARGUS, but it’s easy to implement and performs well across a wide parameter range.
def argus_rvs_numpy(chi: float, size=1, rng: np.random.Generator | None = None) -> np.ndarray:
"""Sample from standardized ARGUS(chi) using only NumPy."""
chi = float(chi)
if chi <= 0:
raise ValueError("chi must be > 0")
if rng is None:
rng = np.random.default_rng()
size1d = tuple(np.atleast_1d(size))
n = int(np.prod(size1d))
out = np.empty(n, dtype=float)
simulated = 0
chi2 = chi * chi
if chi <= 0.5:
# Case 1: propose from g(x) = 3*x*sqrt(1-x^2), accept with exp(-chi^2(1-x^2)/2)
d = -chi2 / 2
while simulated < n:
k = n - simulated
u = rng.uniform(size=k)
v = rng.uniform(size=k)
z = v ** (2 / 3) # z = 1 - x^2 under the proposal
accept = np.log(u) <= d * z
num_accept = int(np.sum(accept))
if num_accept:
out[simulated : simulated + num_accept] = np.sqrt(1 - z[accept])
simulated += num_accept
elif chi <= 1.8:
# Case 2: propose from g(x) ∝ x*exp(-chi^2(1-x^2)/2) and accept with sqrt(1-x^2)
echi = np.exp(-chi2 / 2)
while simulated < n:
k = n - simulated
u = rng.uniform(size=k)
v = rng.uniform(size=k)
# z <= 0, and x = sqrt(1 + z)
z = 2 * np.log(echi * (1 - v) + v) / chi2
accept = (u * u + z) <= 0
num_accept = int(np.sum(accept))
if num_accept:
out[simulated : simulated + num_accept] = np.sqrt(1 + z[accept])
simulated += num_accept
else:
# Case 3: conditional Gamma(3/2, 1) for Y in [0, chi^2/2]
y = np.empty(n, dtype=float)
while simulated < n:
k = n - simulated
g = rng.standard_gamma(shape=1.5, size=k) # Gamma(k=3/2, theta=1)
accept = g <= chi2 / 2
num_accept = int(np.sum(accept))
if num_accept:
y[simulated : simulated + num_accept] = g[accept]
simulated += num_accept
out = np.sqrt(1 - 2 * y / chi2)
return out.reshape(size1d)
# Sanity check: NumPy sampler vs SciPy moments
chi_check = 2.5
s = argus_rvs_numpy(chi_check, size=200_000, rng=rng)
print("sample mean:", s.mean())
print("theory mean:", argus.mean(chi_check))
print("sample var:", s.var())
print("theory var:", argus.var(chi_check))
sample mean: 0.7621546176726476
theory mean: 0.7627785650581257
sample var: 0.035698472048319115
theory var: 0.035554890425441354
8) Visualization (PDF, CDF, Monte Carlo)#
We’ll visualize how (\chi) changes the shape and verify the sampler by overlaying a Monte Carlo histogram on the theoretical PDF.
x = np.linspace(0, 1, 800)
chis_plot = [0.2, 0.7, 2.5, 6.0]
fig = make_subplots(
rows=1,
cols=2,
subplot_titles=("PDF", "CDF"),
)
for chi in chis_plot:
fig.add_trace(go.Scatter(x=x, y=argus_pdf(x, chi), name=f"chi={chi}", mode="lines"), row=1, col=1)
fig.add_trace(go.Scatter(x=x, y=argus_cdf(x, chi), name=f"chi={chi}", mode="lines", showlegend=False), row=1, col=2)
fig.update_xaxes(title_text="x (standardized)", row=1, col=1)
fig.update_xaxes(title_text="x (standardized)", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="CDF", row=1, col=2)
fig.update_layout(title="ARGUS distribution for different chi", width=950)
fig.show()
# Monte Carlo overlay
chi_mc = 2.5
n_mc = 40_000
samples = argus_rvs_numpy(chi_mc, size=n_mc, rng=rng)
hist = np.histogram(samples, bins=70, range=(0, 1), density=True)
bins = hist[1]
centers = 0.5 * (bins[:-1] + bins[1:])
fig2 = go.Figure()
fig2.add_trace(go.Bar(x=centers, y=hist[0], name="MC histogram", opacity=0.4))
fig2.add_trace(go.Scatter(x=x, y=argus_pdf(x, chi_mc), name="theoretical pdf", mode="lines"))
fig2.update_layout(
title=f"Monte Carlo check (chi={chi_mc})",
xaxis_title="x",
yaxis_title="density",
bargap=0.02,
)
fig2.show()
9) SciPy integration (scipy.stats.argus)#
SciPy exposes ARGUS as a standard rv_continuous distribution:
argus.pdf(x, chi)/argus.logpdf(x, chi)argus.cdf(x, chi)/argus.sf(x, chi)argus.rvs(chi, size=..., random_state=...)argus.fit(data, ...)
Because the support is bounded, it’s common to pre-normalize data to ((0,1)) and fit only (\chi) by fixing loc=0, scale=1.
chi = 2.5
x = np.linspace(0, 1, 6)
print("x:", x)
print("pdf:", argus.pdf(x, chi))
print("cdf:", argus.cdf(x, chi))
# Random variates
r = argus.rvs(chi, size=5, random_state=rng)
print("\nrvs:", r)
# Fit chi (fix loc/scale)
chi_true = 2.5
x_data = argus.rvs(chi_true, size=2000, random_state=rng)
chi_hat, loc_hat, scale_hat = argus.fit(x_data, floc=0.0, fscale=1.0)
print("\ntrue chi:", chi_true)
print("fit chi :", chi_hat)
print("(loc, scale fixed to)", loc_hat, scale_hat)
x: [0. 0.2 0.4 0.6 0.8 1. ]
pdf: [0. 0.13515406 0.36789472 0.89991027 2.15877251 0. ]
cdf: [0. 0.01283353 0.06035862 0.17934912 0.46903877 1. ]
rvs: [0.44379125 0.96750725 0.75501664 0.64295312 0.792296 ]
true chi: 2.5
fit chi : 2.528417968750004
(loc, scale fixed to) 0.0 1.0
10) Statistical use cases#
10.1 Hypothesis testing#
A common workflow is to compare a fixed-shape ARGUS background (e.g. (\chi=\chi_0)) against a fitted (\chi) using a likelihood ratio test.
Caveat: the usual (\chi^2) calibration is asymptotic and can be inaccurate for small samples.
10.2 Bayesian modeling#
Treat (\chi) as an unknown parameter with a prior (e.g. log-normal). Because it’s 1D, a grid posterior is often sufficient.
10.3 Generative modeling#
ARGUS is frequently used as a background component inside a mixture model:
[ p(x)=\pi,p_{\text{signal}}(x) + (1-\pi),p_{\text{bkg}}(x),\qquad p_{\text{bkg}}(x)=\text{ARGUS}(\chi). ]
Below is a lightweight demo of all three.
# --- 10.1 Likelihood ratio test demo ---
chi_true = 2.5
x = argus.rvs(chi_true, size=600, random_state=rng)
res = minimize_scalar(lambda c: -argus_loglik(c, x), bounds=(1e-6, 30.0), method="bounded")
chi_hat = float(res.x)
chi0 = 1.0
lr_stat = 2 * (argus_loglik(chi_hat, x) - argus_loglik(chi0, x))
print("chi_hat:", chi_hat)
print("LRT stat:", lr_stat)
print("approx p-value:", chi2.sf(lr_stat, df=1))
# --- 10.2 Bayesian grid posterior for chi ---
# Prior: log chi ~ Normal(mu=0, sigma=1) (=> chi is log-normal)
chis = np.linspace(0.05, 10.0, 400)
log_prior = norm.logpdf(np.log(chis), loc=0.0, scale=1.0) - np.log(chis) # Jacobian for chi -> log chi
log_like = np.array([argus_loglik(c, x) for c in chis])
log_post_unnorm = log_like + log_prior
log_post_unnorm -= np.max(log_post_unnorm)
post = np.exp(log_post_unnorm)
post /= np.trapz(post, chis)
post_cdf = np.cumsum(post)
post_cdf /= post_cdf[-1]
def quantile(q):
return float(np.interp(q, post_cdf, chis))
ci_low, ci_high = quantile(0.05), quantile(0.95)
post_mean = float(np.trapz(chis * post, chis))
print("\nPosterior mean:", post_mean)
print("90% credible interval:", (ci_low, ci_high))
fig = go.Figure(go.Scatter(x=chis, y=post, mode="lines", name="posterior"))
fig.add_vline(x=chi_true, line_dash="dash", line_color="black", annotation_text="true")
fig.add_vrect(x0=ci_low, x1=ci_high, fillcolor="lightblue", opacity=0.3, line_width=0)
fig.update_layout(title="Posterior over chi (grid)", xaxis_title="chi", yaxis_title="density")
fig.show()
# --- 10.3 Simple generative mixture (signal + ARGUS background) ---
# Background: ARGUS
chi_bkg = 6.0
n = 50_000
# Signal: a narrow truncated normal near the endpoint
mu_sig, sigma_sig = 0.97, 0.02
pi_sig = 0.08
n_sig = int(round(pi_sig * n))
n_bkg = n - n_sig
x_bkg = argus_rvs_numpy(chi_bkg, size=n_bkg, rng=rng)
# Truncated normal via rejection (fine for a demo)
x_sig = []
while len(x_sig) < n_sig:
z = rng.normal(loc=mu_sig, scale=sigma_sig, size=n_sig)
z = z[(z > 0) & (z < 1)]
x_sig.extend(z.tolist())
x_sig = np.array(x_sig[:n_sig])
x_mix = np.concatenate([x_bkg, x_sig])
fig = px.histogram(
x_mix,
nbins=120,
histnorm="probability density",
title="Mixture example: ARGUS background + truncated-normal signal",
labels={"value": "x"},
)
xx = np.linspace(0, 1, 800)
fig.add_trace(
go.Scatter(
x=xx,
y=(1 - pi_sig) * argus_pdf(xx, chi_bkg),
name="(1-π) * ARGUS pdf",
mode="lines",
)
)
fig.update_layout(showlegend=True)
fig.show()
chi_hat: 2.526546565812786
LRT stat: 294.2642747018417
approx p-value: 5.853726097616042e-66
Posterior mean: 2.5215534518729363
90% credible interval: (2.3998870486119035, 2.6156776669782578)
11) Pitfalls#
Parameter constraints: (\chi\le 0) is invalid.
Boundary values: the theoretical support is ((0,1)). Real data may contain exact 0/1 due to rounding; for likelihood work you may need to clip slightly, e.g.
x = np.clip(x, 1e-12, 1-1e-12).Numerical stability near 1: use
log1p(-x*x)rather thanlog(1-x**2).Sampling efficiency: the simple truncated-Gamma method is great for large (\chi) but can reject a lot when (\chi) is small; the piecewise sampler avoids that.
Fitting (
fit) surprises: by default, SciPy also fitslocandscale. If your data are already standardized, fix them withfloc=0, fscale=1.
12) Summary#
argusis a continuous distribution on ((0,1)) with shape parameter (\chi>0).Its PDF combines a phase-space term (x\sqrt{1-x^2}) and an exponential tilt controlled by (\chi).
The CDF is especially clean in terms of (\Psi(\chi)): (F(x)=1-\Psi(\chi\sqrt{1-x^2})/\Psi(\chi)).
Mean and second moment have closed forms (Bessel/incomplete-gamma), while skewness/kurtosis are usually computed numerically.
Sampling can be implemented with rejection methods and a truncated-Gamma transformation; SciPy provides a fast and robust implementation in
scipy.stats.argus.